!New version for second article
module fonctions
implicit none
!------------------------------------------------------------------------
double precision, parameter :: x = 0.5d0 !Mesh refinement: x smaller => mesh smaller without syst.size modification
double precision, parameter :: y = 1.0d0 !System size: y bigger => system bigger without mesh modification
!Adjust spatial mesh to prevent instabilities (fine enough), reasonable CPU time (not too fine), and boundaries
!Euler stability criterion: (DX)^2/(D*(Dt)) > 2, D=10^{-9} [m^2/s], DX=realspace, Dt=realtime
!------------------------------------------------------------------------
double precision, parameter :: U=-380.65d0*0.21d0!Imposed potential: Nl*y => U*y
double precision, parameter :: deltat = 0.0005d0 !time discretisation.
integer, parameter :: Nsavedata = 960000 !Saving data every Nsavedata iterations
integer, parameter :: Nsavedata2 = 96 !Saving NEW data every Nsavedata it. !Not more than 99 saves, or last one in "liese-R10.txt" overwritten
integer, parameter :: nsave = 15 !maximum number saves for Nsavedata. Then last one overwritten. Alonissos: =2. Amorgos: =15.
!------------------------------------------------------------------------
double precision, parameter :: realspace = .01d0*0.5d0 !Conversion factor for real space [centimeter]: /x
double precision, parameter :: realtime = (10.d0/60.d0)*(0.5d0)**2 !Conversion factor for real time [minutes]: /x^2
integer, parameter :: Nl = ((100.d0/0.5d0)*1.0d0) !#sites to left; mesh size: Nl*realspace=constant | /x
integer, parameter :: Nr = ((2000.d0/0.5d0)*1.0d0) !#sites to right; mesh size: keep Nr*realspace=constant | /x
!------------------------------------------------------------------------
double precision, parameter :: a0 = 10.d0 !Initial concentration of A
double precision, parameter :: b0 = 0.1d0 !Initial concentration of B
double precision, dimension(1:4), parameter :: D = (/ 1.0d0, 1.d0, 1.d0, 1.d0 /) !Diffusion constants for A+,A-,B+,B-
!------------------------------------------------------------------------
logical, parameter :: infiniterate = .false. !rate of reaction: infinite or finite
double precision, parameter :: rate = 10.d0 !rate of reaction if finite
!------------------------------------------------------------------------
character*(*), parameter :: filename0 = 'liese-state.txt'
character*(*), parameter :: filename1 = 'liese-C.txt' !concentration of C(x,t): t,{xj,C(xj)}
character*(*), parameter :: filename2 = 'liese-R.txt' !reaction rate R(x,t): t,{xj,R(xj)}
character*(*), parameter :: filename3 = 'liese-I.txt' !current I(t): t,I(t)
character*(*), parameter :: filename4 = 'liese-E.txt' !electric field E(t): t,{xj,E(xj)}
character*(*), parameter :: filename4b = 'liese-Et.txt' !electric field E(t) at the front: t,E(xf)
character*(*), parameter :: filename5 = 'liese-ap.txt' !concentration of A+: t,{xj,A+(xj)}
character*(*), parameter :: filename5b = 'liese-apt.txt' !concentration of A+ at the front: t,A+(xfront)
character*(*), parameter :: filename6 = 'liese-am.txt' !concentration of A-: t,{xj,A-(xj)}
character*(*), parameter :: filename6b = 'liese-amt.txt' !concentration of A- at the front: t,A-(xfront)
character*(*), parameter :: filename7 = 'liese-bp.txt' !concentration of B+: t,{xj,B+(xj)}
character*(*), parameter :: filename7b = 'liese-bpt.txt' !concentration of B+ at the front: t,B+(xfront)
character*(*), parameter :: filename8 = 'liese-bm.txt' !concentration of B-: t,{xj,B-(xj)}
character*(*), parameter :: filename8c = 'liese-bmt.txt' !concentration of B- at the front: t,B-(xfront)
character*(*), parameter :: filename8b = 'liese-AB.txt' !concentration of all species-: t,{xj,all(xj)}
character*(*), parameter :: filename9 = 'liese-q.txt' !charge: t,{xj,(A+ - A- + B+ - B-)(xj)}
character*(*), parameter :: filename10 = 'liese-Qtot.txt' !total charge Q(t): t,Q(t)
character*(*), parameter :: filename11 = 'liese-Cfront.txt' !position of the C front: t,Cfront(t)
character*(*), parameter :: filenamenewa = 'liese-w2wrmaxt.txt' !w^2(t), w(t), R_{max}(t)
character*(*), parameter :: filenamenewb1 = 'liese-R1.txt' !reaction rate R(x,t): t,{xj,R(xj)}, 0-9
character*(*), parameter :: filenamenewb2 = 'liese-R2.txt' !reaction rate R(x,t): t,{xj,R(xj)}, 10-19
character*(*), parameter :: filenamenewb3 = 'liese-R3.txt' !reaction rate R(x,t): t,{xj,R(xj)}, 20-29
character*(*), parameter :: filenamenewb4 = 'liese-R4.txt' !reaction rate R(x,t): t,{xj,R(xj)}, 30-39
character*(*), parameter :: filenamenewb5 = 'liese-R5.txt' !reaction rate R(x,t): t,{xj,R(xj)}, 40-49
character*(*), parameter :: filenamenewb6 = 'liese-R6.txt' !reaction rate R(x,t): t,{xj,R(xj)}, 50-59
character*(*), parameter :: filenamenewb7 = 'liese-R7.txt' !reaction rate R(x,t): t,{xj,R(xj)}, 60-69
character*(*), parameter :: filenamenewb8 = 'liese-R8.txt' !reaction rate R(x,t): t,{xj,R(xj)}, 70-79
character*(*), parameter :: filenamenewb9 = 'liese-R9.txt' !reaction rate R(x,t): t,{xj,R(xj)}, 80-89
character*(*), parameter :: filenamenewb10 = 'liese-R10.txt' !reaction rate R(x,t): t,{xj,R(xj)}, 90-99
double precision, parameter :: Cfront = 0.0001 !If C < Cfront then this is the position of the C front
integer, parameter :: Ns = ((2100.d0/0.5d0)*1.0d0) !total number of sites. Do not adjust.
real, dimension(1:nsave,1:Ns) :: tmpC,tmpE,tmpR,tmpAP,tmpAM,tmpBP,tmpBM,tmpQ,tmpAB !temporary for data saving.
real, dimension(1:10,1:Ns) :: tmpR2
integer :: pouet
contains

!------------------------------------------------------------------------
subroutine DerivativesRHS(AB,E,Dxx,Ex,XI)
double precision, dimension(1:4,1:Ns), intent(in) :: AB
double precision, dimension(1:Ns), intent(inout) :: E
double precision, dimension(1:4,1:Ns), intent(inout) :: Ex,Dxx
double precision, intent(out) :: XI
double precision, dimension(1:Ns) :: sum1,sum2 !Sums in expr. of E
double precision :: tmp
integer :: i
!E field
do i=2,Ns-1
	sum1(i) = D(1)*(AB(1,i+1)-AB(1,i)) - D(2)*(AB(2,i+1)-AB(2,i)) &
					+ D(3)*(AB(3,i+1)-AB(3,i)) - D(4)*(AB(4,i+1)-AB(4,i))
	sum2(i) = (D(1)*(AB(1,i+1)+AB(1,i)) + D(2)*(AB(2,i+1)+AB(2,i)) &
					+ D(3)*(AB(3,i+1)+AB(3,i)) + D(4)*(AB(4,i+1)+AB(4,i)))/2.d0
end do
sum1(1)=0.d0
sum2(1)=a0*(D(1)+D(2))
sum1(Ns)=0.d0
sum2(Ns)=b0*(D(3)+D(4))
tmp=0.d0
do i=1,Ns
	tmp = tmp + 1.d0/sum2(i)
end do
XI=U/tmp
E(:) = (XI+sum1(:))/sum2(:)
!Right-hand side: $d^2 c_i / d x^2$ (Dxx) and $-z_i d / dx c_i E$ (Ex)
do i=2,Ns-1
	Ex(1:4:2,i) = - (E(i)*(AB(1:4:2,i+1)+AB(1:4:2,i)) - E(i-1)*(AB(1:4:2,i)+AB(1:4:2,i-1)))/2.d0
	Ex(2:4:2,i) = + (E(i)*(AB(2:4:2,i+1)+AB(2:4:2,i)) - E(i-1)*(AB(2:4:2,i)+AB(2:4:2,i-1)))/2.d0
	Dxx(:,i) = AB(:,i+1) - 2.d0*AB(:,i) + AB(:,i-1)
end do
return
end subroutine DerivativesRHS
!------------------------------------------------------------------------
subroutine RK(AB,E,Dxx,Ex)
double precision, dimension(1:4,1:Ns), intent(inout) :: AB
double precision, dimension(1:Ns), intent(in) :: E
double precision, dimension(1:4,1:Ns), intent(in) :: Ex,Dxx
double precision, dimension(1:4,1:Ns) :: tmpAB,p1,p2,p3,p4,tmpEx,tmpDxx
double precision, dimension(1:Ns) :: tmpE
double precision :: tmp
integer :: i
	do i=1,4
		p1(i,:) = D(i)*(Dxx(i,:)+Ex(i,:))
	end do
	tmpAB = AB+p1*deltat/2.d0
!	tmpAB(1,:) = tmpAB(2,:) + tmpAB(4,:) - tmpAB(3,:) !a-posteriori correction
	tmpE = E
	tmpEx = Ex
	tmpDxx = Dxx
	call DerivativesRHS(tmpAB,tmpE,tmpDxx,tmpEx,tmp)
	do i=1,4
		p2(i,:) = D(i)*(tmpDxx(i,:)+tmpEx(i,:))
	end do
	tmpAB = AB+p2*deltat/2.d0
!	tmpAB(1,:) = tmpAB(2,:) + tmpAB(4,:) - tmpAB(3,:) !a-posteriori correction
	call DerivativesRHS(tmpAB,tmpE,tmpDxx,tmpEx,tmp)
	do i=1,4
		p3(i,:) = D(i)*(tmpDxx(i,:)+tmpEx(i,:))
	end do
	tmpAB = AB+p3*deltat
!	tmpAB(1,:) = tmpAB(2,:) + tmpAB(4,:) - tmpAB(3,:) !a-posteriori correction
	call DerivativesRHS(tmpAB,tmpE,tmpDxx,tmpEx,tmp)
	do i=1,4
		p4(i,:) = D(i)*(tmpDxx(i,:)+tmpEx(i,:))
	end do
	AB = AB + deltat/6.d0*(p1+2.d0*p2+2.d0*p3+p4)
!	AB(1,:) = AB(2,:) + AB(4,:) - AB(3,:) !a-posteriori correction
!	where (AB < 0.d0) AB = 0.d0
return
end subroutine RK
!------------------------------------------------------------------------
subroutine Euler(AB,Dxx,Ex)
double precision, dimension(1:4,1:Ns),intent(inout) :: AB
double precision, dimension(1:4,1:Ns), intent(in) :: Dxx,Ex
integer :: i
	do i=1,4
		AB(i,:) = AB(i,:) + deltat*D(i)*(Dxx(i,:)+Ex(i,:))
	end do
return
end subroutine Euler
!----------------------------------------------------------------------
subroutine erasefiles
	open(unit=1,file=filename0)
	close(1,status='DELETE')
	open(unit=1,file=filename1)
	close(1,status='DELETE')
	open(unit=1,file=filename2)
	close(1,status='DELETE')
	open(unit=1,file=filename3)
	close(1,status='DELETE')
	open(unit=1,file=filename4)
	close(1,status='DELETE')
	open(unit=1,file=filename4b)
	close(1,status='DELETE')
	open(unit=1,file=filename5)
	close(1,status='DELETE')
	open(unit=1,file=filename5b)
	close(1,status='DELETE')
	open(unit=1,file=filename6)
	close(1,status='DELETE')
	open(unit=1,file=filename6b)
	close(1,status='DELETE')
	open(unit=1,file=filename7)
	close(1,status='DELETE')
	open(unit=1,file=filename7b)
	close(1,status='DELETE')
	open(unit=1,file=filename8)
	close(1,status='DELETE')
	open(unit=1,file=filename8c)
	close(1,status='DELETE')
	open(unit=1,file=filename8b)
	close(1,status='DELETE')
	open(unit=1,file=filename9)
	close(1,status='DELETE')
	open(unit=1,file=filename10)
	close(1,status='DELETE')
	open(unit=1,file=filename11)
	close(1,status='DELETE')
return
end subroutine erasefiles
!----------------------------------------------------------------------
subroutine erasefiles2
	open(unit=1,file=filenamenewa)
	close(1,status='DELETE')
	open(unit=1,file=filenamenewb1)
	close(1,status='DELETE')
	open(unit=1,file=filenamenewb2)
	close(1,status='DELETE')
	open(unit=1,file=filenamenewb3)
	close(1,status='DELETE')
	open(unit=1,file=filenamenewb4)
	close(1,status='DELETE')
	open(unit=1,file=filenamenewb5)
	close(1,status='DELETE')
	open(unit=1,file=filenamenewb6)
	close(1,status='DELETE')
	open(unit=1,file=filenamenewb7)
	close(1,status='DELETE')
	open(unit=1,file=filenamenewb8)
	close(1,status='DELETE')
	open(unit=1,file=filenamenewb9)
	close(1,status='DELETE')
	open(unit=1,file=filenamenewb10)
	close(1,status='DELETE')
return
end subroutine erasefiles2
!----------------------------------------------------------------------
subroutine savedata(C,R,XI,E,AB,k,time,front)
double precision, dimension(1:Ns), intent(in) :: C,R,E
double precision, intent(in) :: XI
double precision, dimension(1:4,1:Ns), intent(in) :: AB
integer, intent(in) :: k
double precision, intent(in) :: time
real :: Nlrsp
integer :: i,j,n,front
double precision :: minR
	open(unit=1,file=filename0,position='APPEND')
	write(1,'("Number of save: ",i5,"     Iteration: ",i12,"    Time [min]: ",f12.2)') k,(k-1)*Nsavedata,time*realtime
	close(1)
	n=k
	if (n >= nsave) n=nsave !if more than nsave steps saved, data number nsave overwritten
	!temporary saving
	tmpC(n,:) = real(C(:))
	tmpR(n,:) = real(R(:)/deltat)
	tmpE(n,:) = real(E(:))
	tmpAP(n,:) = real(AB(1,:))
	tmpAM(n,:) = real(AB(2,:))
	tmpBP(n,:) = real(AB(3,:))
	tmpBM(n,:) = real(AB(4,:))
	tmpQ(n,:) = real(AB(1,:)-AB(2,:)+AB(3,:)-AB(4,:))
	tmpAB(n,:) = real(AB(1,:)+AB(2,:)+AB(3,:)+AB(4,:))
	!saving in columns for GnuPlot
	Nlrsp = real(Nl)*real(realspace)
	
	!position of the front defined by the maximum of the production rate R
	minR = 0.d0
	do i=1,Nr
		if (R(i+Nl) > minR) then
			minR = R(i+Nl)
			front = i
		end if
	end do	
	
	open(unit=1,file=filename11,position='APPEND')
	write(1,*) real(time*realtime),real(dble(front)*realspace)
	close(1) !the position of the front (in cm) as a function of time (in
	! min.)
	
	open(unit=1,file=filename1)
	do i=1,Ns
		write(1,*) (-Nlrsp+real(dble(i)*realspace),tmpC(j,i),j=1,n)
	end do !concentration in molar
	close(1)
	
	open(unit=1,file=filename2)
	do i=1,Ns
		write(1,*) (-Nlrsp+real(dble(i)*realspace),tmpR(j,i),j=1,n)
	end do
	close(1)

	open(unit=1,file=filename3,position='APPEND')
	write(1,*) real(time*realtime),XI*realspace/(realtime*60.d0)*96.37d0
	close(1) !current in A/cm^2
	
	open(unit=1,file=filename4)
	do i=1,Ns
	write(1,*)(-Nlrsp+real(dble(i)*realspace),tmpE(j,i)/(38.65d0*realspace*0.01d0),j=1,n)
	end do ! local field in V/m
	
	close(1)
	open(unit=1,file=filename4b,position='APPEND')
	write(1,*)real(time*realtime),minval(E(:))/(38.65d0*realspace*0.01d0)
	close(1) !minvalue of the field in V/m
	
	open(unit=1,file=filename5)
	do i=1,Ns
		write(1,*) (-Nlrsp+real(dble(i)*realspace),tmpAP(j,i),j=1,n)
	end do
	close(1)
	
	open(unit=1,file=filename5b,position='APPEND')
	write(1,*)real(time*realtime),AB(1,Nl+front)
	close(1)
	
	open(unit=1,file=filename6)
	do i=1,Ns
		write(1,*) (-Nlrsp+real(dble(i)*realspace),tmpAM(j,i),j=1,n)
	end do
	close(1)
	
	open(unit=1,file=filename6b,position='APPEND')
	write(1,*) real(time*realtime),AB(2,Nl+front)
	close(1)
	
	open(unit=1,file=filename7)
	do i=1,Ns
		write(1,*) (-Nlrsp+real(dble(i)*realspace),tmpBP(j,i),j=1,n)
	end do
	close(1)
	
	open(unit=1,file=filename7b,position='APPEND')
	write(1,*) real(time*realtime),AB(3,Nl+front)
	close(1)
	
	open(unit=1,file=filename8)
	do i=1,Ns
		write(1,*) (-Nlrsp+real(dble(i)*realspace),tmpBM(j,i),j=1,n)
	end do
	close(1)
	
	open(unit=1,file=filename8c,position='APPEND')
	write(1,*) real(time*realtime),AB(4,Nl+front)
	close(1)
	
	open(unit=1,file=filename8b)
	do i=1,Ns
		write(1,*) (-Nlrsp+real(dble(i)*realspace),tmpAB(j,i),j=1,n)
	end do
	close(1)
	
	open(unit=1,file=filename9)
	do i=1,Ns
		write(1,*)(-Nlrsp+real(dble(i)*realspace),tmpQ(j,i)/(realspace)**2,j=1,n)
	end do
	close(1)
	
	open(unit=1,file=filename10,position='APPEND')
	write(1,*)real(time*realtime),(sum(AB(1,:))-sum(AB(2,:))+sum(AB(3,:))-sum(AB(4,:)))/(realspace)**2
	close(1)
return
end subroutine savedata

!----------------------------------------------------------------------
subroutine savedata2(R,k2,time)
double precision, dimension(1:Ns), intent(in) :: R
double precision, intent(in) :: time
integer, intent(in) :: k2
real :: Nlrsp
integer :: i,j,front,n,m
double precision :: minR,w2,den
	n=mod(k2-1,10) !unités
	m=(k2-1-n)/10 !dixaines

	!temporary saving
	if (m<10) tmpR2(n+1,:) = real(R(:)/deltat)
	if (m>9) tmpR2(10,:) = real(R(:)/deltat)
	!saving in columns for GnuPlot
	Nlrsp = real(Nl)*real(realspace)
	
	!position of the front defined by the maximum of the production rate R
	minR = 0.d0
	do i=1,Nr
		if (R(i+Nl) > minR) then
			minR = R(i+Nl)
			front = i+Nl
		end if
	end do

	!Calculation of the integral
  den=0.d0
	do i=1,Ns-1
		den=den + R(i)
	end do
  w2=0.d0
  do i=1,Ns-1
    w2=w2+dble((i-front)*real(realspace))**2*R(i)
  end do
  w2=w2/den

	open(unit=1,file=filenamenewa,position='APPEND')
	write(1,*) real(time*realtime),w2,dsqrt(w2),minR
	close(1)
	
	if (m==0) then
		open(unit=1,file=filenamenewb1)
		do i=1,Ns
			write(1,*) (-Nlrsp+real(dble(i)*realspace),tmpR2(j,i),j=1,n+1)
		end do
		close(1)
!  	write(*,*) "Save nb: ",k2," File: ",m+1," Position in file: ",n+1
	end if
	if (m==1) then
		open(unit=1,file=filenamenewb2)
		do i=1,Ns
			write(1,*) (-Nlrsp+real(dble(i)*realspace),tmpR2(j,i),j=1,n+1)
		end do
		close(1)
	end if
	if (m==2) then
		open(unit=1,file=filenamenewb3)
		do i=1,Ns
			write(1,*) (-Nlrsp+real(dble(i)*realspace),tmpR2(j,i),j=1,n+1)
		end do
		close(1)
	end if
	if (m==3) then
		open(unit=1,file=filenamenewb4)
		do i=1,Ns
			write(1,*) (-Nlrsp+real(dble(i)*realspace),tmpR2(j,i),j=1,n+1)
		end do
		close(1)
	end if
	if (m==4) then
		open(unit=1,file=filenamenewb5)
		do i=1,Ns
			write(1,*) (-Nlrsp+real(dble(i)*realspace),tmpR2(j,i),j=1,n+1)
		end do
		close(1)
	end if	
	if (m==5) then
		open(unit=1,file=filenamenewb6)
		do i=1,Ns
			write(1,*) (-Nlrsp+real(dble(i)*realspace),tmpR2(j,i),j=1,n+1)
		end do
		close(1)
	end if	
	if (m==6) then
		open(unit=1,file=filenamenewb7)
		do i=1,Ns
			write(1,*) (-Nlrsp+real(dble(i)*realspace),tmpR2(j,i),j=1,n+1)
		end do
		close(1)
	end if	
	if (m==7) then
		open(unit=1,file=filenamenewb8)
		do i=1,Ns
			write(1,*) (-Nlrsp+real(dble(i)*realspace),tmpR2(j,i),j=1,n+1)
		end do
		close(1)
	end if	
	if (m==8) then
		open(unit=1,file=filenamenewb9)
		do i=1,Ns
			write(1,*) (-Nlrsp+real(dble(i)*realspace),tmpR2(j,i),j=1,n+1)
		end do
		close(1)
	end if
	if (m==9) then
		open(unit=1,file=filenamenewb10)
		do i=1,Ns
			write(1,*) (-Nlrsp+real(dble(i)*realspace),tmpR2(j,i),j=1,n+1)
		end do
		close(1)
	end if
	if (m>9) then
		open(unit=1,file=filenamenewb10)
		do i=1,Ns
			write(1,*) (-Nlrsp+real(dble(i)*realspace),tmpR2(j,i),j=1,10)
		end do
		close(1)
	end if
	
return
end subroutine savedata2
!----------------------------------------------------------------------
end module fonctions

!----------------------------------------------------------------------
!----------------------------------------------------------------------

program liese
use fonctions
implicit none
double precision, dimension(1:4,1:Ns) :: AB !A+,A-,B+,B-,C
double precision, dimension(1:4,1:Ns) :: Dxx !Space derivatives of A+,A-,B+,B-
double precision, dimension(1:4,1:Ns) :: Ex !Derivatives of term with elect. field for A+,A-,B+,B-
double precision, dimension(1:Ns) :: E,R,C !electrical field E, reaction term, C
double precision :: time,XI !time t, current I
integer :: i,j,k,front,i2,k2,n
!Initial conditions: square profile
AB(1:2,1:Nl)=a0
AB(3:4,1:Nl)=0.d0
AB(1:2,Nl+1:Ns)=0.d0
AB(3:4,Nl+1:Ns)=b0
C=0.d0
R=0.d0
Ex=0.d0
Dxx=0.d0
call DerivativesRHS(AB,E,Ex,Dxx,XI)
time=0.d0
i=0
i2=0
k=1
k2=1
pouet=0
call erasefiles
call erasefiles2
mainloop: do
	i=i+1
	i2=i2+1
	time=time+deltat
!	call Euler(AB,Dxx,Ex)
	call RK(AB,E,Dxx,Ex)
		do j=1,Ns
!		if (infiniterate .eqv. .true.) then
!			R(j) = min(AB(2,j),AB(3,j)) !for infinite reaction rate
!		else
!			if (i==1 .and. k==1) R(j) = rate*deltat*AB(2,j)*AB(3,j)
!			if (R(j) > min(AB(2,j),AB(3,j))) R(j) = min(AB(2,j),AB(3,j)) !for finite reaction rate
!		end if
!		AB(2,j) = AB(2,j) - R(j)
!		AB(3,j) = AB(3,j) - R(j)
		if (infiniterate .eqv. .false.) then
		R(j) = rate*deltat*AB(2,j)*AB(3,j) !for finite reaction rate
		if (R(j) > min(AB(2,j),AB(3,j))) R(j) = min(AB(2,j),AB(3,j)) !for finite reaction rate
	        end if
		AB(2,j) = AB(2,j) - R(j)
		AB(3,j) = AB(3,j) - R(j)
		
	end do
	
	C = C + R
	call DerivativesRHS(AB,E,Dxx,Ex,XI)
	if (i==Nsavedata .or. (i==1 .and. k==1)) then
		call savedata(C,R,XI,E,AB,k,time,front)
		if (k /= 1) i=0
		k=k+1
	end if
	if (i2==Nsavedata2 .or. (i2==1 .and. k2==1)) then
		call savedata2(R,k2,time)
		if (k2 /= 1) i2=0
		k2=k2+1
	end if	
	
	if (time >= 1800.0d0/(0.5d0)**2) exit mainloop
end do mainloop
end program liese
